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Abstract 

The conditions that can lead to the exploitative depletion of a shared resource, i.e, the tragedy 
of the commons, can be reformulated as a game of prisoner's dilemma: while preserving the common 
resource is in the best interest of the group, over-consumption is in the interest of each particular 
individual at any given point in time. One way to try and prevent the tragedy of the commons 
is through infliction of punishment for over-consumption and/or encouraging under-consumption, 
thus selecting against over-consumers. Here, the effectiveness of various punishment functions in an 
evolving consumer-resource system is evaluated within a framework of a parametrically heterogeneous 
system of ordinary differential equations (ODEs). Conditions leading to the possibility of sustainable 
coexistence with the common resource for a subset of cases are identified analytically using adaptive 
dynamics; the effects of punishment on heterogeneous populations with different initial composition 
are evaluated using the Reduction theorem for replicator equations. Obtained results suggest that 
one cannot prevent the tragedy of the commons through rewarding of under-consumers alone - there 
must also be an implementation of some degree of punishment that increases in a non-linear fashion 
with respect to over-consumption and which may vary depending on the initial distribution of clones 
in the population. 

Keywords: tragedy of the commons, mathematical model, Reduction theorem, adaptive dynam- 
ics 
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Introduction 



Ecological and social systems are complex and adaptive, composed of diverse agents that are intercon- 
nected and interdependent |121 1151 116| . Heterogeneity within these systems often drives the evolution 
and adaptability of the system components enabling them to withstand and recover from environmental 
perturbations. However, it is also heterogeneity that makes the appearance and short-term prosper- 
ity of over-consumers possible, which in turn can eventually lead to exhaustion of the common resources 
(tragedy of the commons |4J) and consequent collapse of the entire population, also known as evolutionary 
suicide [3]). 

Elinor Ostrom has focused on the question of avoiding the tragedy of the commons from the point of 
view of collective decision making in small fisheries [2] . She observed that the mutually satisfactory and 
functioning institution of collective action could be developed in small communities where the possible 
resource over-consumption by each individual was immediately noticeable and punished, such as in water- 
monitoring systems instituted in Spanish huertas, where each user was able to closely watch the other, 
and where punishment for stealing scarce water was instituted not only in monetary fines but also 
in public humiliation. This illustrates a first path to prevent the tragedy of the commons: infliction of 
punishment tax for over-consumption, effectively selecting against over-consumers. For instance, Ostrom 
cites as an example of a successful adaptive governance system rural villages in Japan, where about 
3 million hectares of forests and mountain meadows were successfully managed even until now in no 
small part because "accounts were kept about who contributed to what to make sure no household 
evaded its responsibilities unnoticed. Only illness, family tragedy, or the absence of able-bodied adults 
... were recognized as excuses for getting out of collective labor... But if there was no acceptable excuse, 
punishment was in order" . Another situation, when the tragedy was successfully avoided, is when the 
community introduces some kind of "social currency", where an individual is rewarded for cooperation 
with social status [TT| ITS] . This is an example of a second approach to preventing the tragedy of the 
commons: bestowing reward/subsidy to under-consumers. 

The dynamics of trade-offs between personal and population good have been studied through classical 
game theory |13 | fT7 1 l5]. However, many of the potentially relevant results have been obtained for largely 
homogeneous populations of players. In this paper we approach the question of preventing the tragedy 
of the commons through instituting punishment/reward systems in heterogeneous populations using two 
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recently developed methods for studying evolving heterogeneous populations in equation-based models, 
namely adaptive dynamics [1] and the Reduction theorem for replicator equations [8] . 

This paper is organized as follows: first we give a derivation of a previously studied model, modified to 
incorporate the effects of a punishment /reward system on overall population survival. Next, we describe 
different approaches for studying the effects of punishment and reward on population composition, namely, 
adaptive dynamics and the Reduction theorem for replicator equations. Then we apply the two techniques 
to a variety of punishment /reward structures, obtaining both analytical and numerical results for sample 
populations that differ in their initial composition with respect to over- and under- consumers. We end 
the paper with a reformulation of the obtained results in the context of the prisoner's dilemma, as well 
as with a comparison of the two modeling methods. 

Model Description 

In order to answer the questions posed in the introduction, we first turn to a modified version of a math- 
ematical model of 'niche construction' proposed by Krakauer et al.|10|. where a population of consumers 
x(t) interact with a collectively constructed 'niche', represented as a common dynamical carrying capacity 
z(t). Rather than restricting the model to a specific physical resource, we take the 'niche' to be a gener- 
alized abstract 'resource', such as available nutrients, energy and other characteristics of the environment 
that may affect both the individuals' survival and fecundity (in case of financial systems, the 'resource' 
could be an inter-convertible unit such as money) that each individual can invest or subtract from the 
'common good'. Different individuals can contribute to the common dynamical carrying capacity, or they 
can subtract from it at different rates. Since sufficient reduction of the common 'carrying capacity' can 
indeed lead to population collapse, this model allows to model effectively in a conceptual framework the 
question of the effects of over-consumption on the survival of the population, 

Within the frameworks of the model, each individual is characterized by his or her own intrinsic value 
of resource consumption, denoted by parameter c; a set of consumers that are characterized by the same 
value of c is referred to as c-clone. The total population of all consumers is N(t) = x c if c takes on 
discrete values and N(t) = j A x c dc if c is continuous, where A denotes the range of possible values of 
parameter c. 

The individuals x c (t) grow according to a logistic growth function, where the carrying capacity is not 
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constant but is a dynamic variable. Each individual invests in resource restoration at a rate proportional 
to (1 — c); an individual is considered to be an over-consumer if c > 1. Moreover, each individual is 
punished or rewarded according to a general continuous function /(c) <E C 1 (K + ) that directly affects the 
fitness of each consumer-producer such that each clone is immediately punished for over-consumption if 
c > 1 or rewarded for under-consumption when c < 1. 

The dynamics of the resource is determined by a natural restoration rate 7 and decay /loss rate dz(t). 
Restoration process is accounted for by the term e(l — c) p^r^TI ) wnere z{t) is restored both proportionally 
to each individual's investment, accounted for by parameter c, and to the total amount of 'resource' 
consumed by the entire population, represented by N(t)/ z(t). 

This yields the following system of equations: 



clones 
resource 



rx c (t) 

intrinsic growth 

7 + 

restoration 



N(t) 



consumption 



kz(t) 

dynamical carrying capacity/ 

EaScW(I-c) 



^W/M , (i) 

punishment / reward 



z(t) + N(t) 



Sz(t) 

decay /inflation 



change in resource caused by consumers 

The meaning of all the variables and parameters summarized in Table [T] 





Meaning 


Range 


z{t) 


abstract generalized 'resource' 


z > 


c 


rate of resource consumption 


c G [c ,Cf] 


x c (t) 


population of individuals competing for the resource 


x c >0 


N(t) 


total population size N(t) — J A x c (t)dc 


N(t) > 


r 


individual proliferation rate 


r > 


k 


'resource' conversion factor 


k > 


e 


efficiency of construction of common 'niche' 


e > 


7 


intrinsic rate of 'resource' restoration independent of x c 


7 > 


S 


per capita rate of natural 'resource' decay 


s > 


Co 


lower boundary value of parameter c 


c > 


c f 


upper boundary value of parameter c 


c f >0 



Table 1: Summary of variables and parameters used in System [T] 
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Table 2: Sample parameter values. 



The case, where f(c) — was previously completely studied in jS], both for the case of a para- 
metrically homogeneous and parametrically heterogeneous system with respect to parameter c. Several 
transitional regimes were identified with increasing parameter c, ranging from sustainable coexistence 
with the common resource with ever decreasing domain of attraction to sustained oscillatory regimes to 
collapse due to exhaustion of the common resource. The results are summarized in Figure [T] In Domain 
1, when the parameter of resource consumption is small, the common carrying capacity remains large, 
successfully supporting the entire population, since no individual is taking more resource than they re- 
store. In Domain 2, a parabolic sector appears near the origin, decreasing the domain of attraction of the 
non-trivial equilibrium point A. The population can still sustainably coexist with the resource even with 
moderate levels of over-consumption but the range of initial conditions, where it is possible, decreases. 
As the value of c is further increased, the range of possible parameter values that allow sustainable co- 
existence with the common resource decreases and is now bounded by the unstable limit cycle, which 
appears around point A through a catastrophic Hopf bifurcation in Domain 3, and via generalized Hopf 
bifurcation in Domain 6. Finally, in Domain 4 and 5, population extinction is inevitable due to extremely 
high over-consumption rates unsupportable by the resource. 

In this paper we will investigate the question of whether timely infliction of punishment for over- 
consumption or rewarding under-consumption can prevent such a collapse and consequently prevent the 
tragedy of the commons. We will also investigate the degree of effectiveness of punishment depending on 
initial composition of the population with respect to parameter c. Finally, we will try to investigate the 
question of whether punishing those who over-consume or rewarding those who under-consume will yield 
better results. These questions will be addressed using two recently developed methods for modeling 
parametrically heterogeneous populations, namely, adaptive dynamics [1 and Reduction theorem for 
replicator equations [5]. 
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Figure 1: Bifurcation diagram of System (JTJ when /(c) = in the (7,0) and (N,z) phase spaces for 
fixed positive parameters e and 5. The non-trivial equilibrium point A is globally stable in Domain 1; it 
shares basins of attraction with equilibrium O in Domains 2 and 3. The separatrix of O and the unstable 
limit cycle that contains point A, serve, correspondingly, as the boundaries of the basins of attraction. 
Only equilibrium O is globally stable in Domains 4, which also contains also unstable non-trivial A, 
and 5, which contains the elliptic sector. Domain 6 exists only for S > 5 + x/24, where the stable limit 
cycle that is in turn contained inside an unstable limit cycle, shares basins of attraction with equilibrium 
O. Boundaries between Domains K, S, H, Nul, C correspond respectively to appearance of an attractive 
sector in a neighborhood of O, appearance of unstable limit cycle containing A, change of stability of 
equilibrium A via Hopf bifurcations, disappearance of positive A and saddle-node bifurcation of limit 
cycles. 
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Adaptive dynamics 

Adaptive dynamics is a series of techniques that have been recently developed to address questions of 
system invasibility by rare 'mutant' clones. The main focus of this method is in evaluating a "mutant's" 
ability to proliferate in an environment set by the 'resident' population pQ. We use this method to evaluate 
different punishment /reward functions and their efficiency in preventing the tragedy of the commons. It is 
worth noting that a more general theory of systems with inheritance, which allows studying the problem of 
'external stability' of a population of 'residents' to invasion by a 'mutant' population, has been developed 
in |2|; it was later published in English in [3]. 

Consider the equation for the dynamics of a rare mutant x m in an environment set by the resident 
x res . The total population size is x m + x res sa x res , since x m is assumed to be present at such low 
frequency that its contribution to the size of the entire population is negligible. 

Let x* es satisfy da J" — , which implies that x* es = ^r(/(c res ) + rc res ). Now introduce a mutant, 
such that 



dx r 



*£m %res • 



(3) 



dt ' kz(t) 

When the two subpopulations interact, the outcome of their interaction is determined by the sign of 
the dominant eigenvalue of System 



dt 



es V ^res 



kz(t) 



dt 
d?. 
dt 



TXm (c< 

7 + e 



A " m I ^ res ) _ 1— x fir \ 

kz{t) > +Xm ^ Cm > 



(4) 



5z 



The stability of the single positive non-trivial equilibrium point (x*, x* m , z*) under the assumption 
that the mutant is rare, i.e., when x* m 0, is determined by the dominant eigenvalue of System Q, 
which is given by A — r(c m — c res ) — f{c res ) + /(c m ): the mutant will be able to invade if A > and it 
will not succeed if A < 0. Invasion fitness of the mutant, i.e., the expected growth rate of a mutant in an 
environment set by the resident, is given by 



ivm 

T-s-ooT 



kz{t) 



f(Cn 



dt = r{c n 



1 ,x 



( )) + /(c m ) = r(c m , E res 



), 



(5) 
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where (^f^) = lim ^ J Q T x, k " z s ( j t ) > dt. The selection gradient, which is defined as the slope of invasioga- 
men fitness, determines if the invasion will be successful (positive sign of the selection gradient predicts 
successful invasion) and is then given by 



D{c m ) = - — r(c m ,£7 res )|[ Cm=Cre j = r + f'(c m ). (6) 

These conditions allow answering a question of "long-term" invasibility, i.e., whether a mutant with a 
slightly difference value of c m will be able to permanently invade the population of individuals, character- 



ized by parameter c res . The points where selection gradient is zero are known as evolutionarily singular 
strategies (ESS) and are denoted here as c* es . Stability of c* es for different types of punishment /reward 
functions is discussed below; summary of some of the possible types of c* es is given in Table [3] 





Adaptive dynamics 


Reduction theorem 


Goal 


Model evolution of parametrically heterogeneous populations 


Assumptions 


Clonal reproduction 




Separation of evolutionary and ecological time scales 


Environment 


Can be variable and affected by changes in 
seascape, or a dancing landscape) 


population composition (also known as 


Population 


Two types: invader and resident 


Any number of types 




Small initial frequency of the invader 


Initial population composition can be 
arbitrary 




Can introduce new "mutants" 


All types must initially be present in the 
population, even if at a near-zero 
frequency 




Requires "Lotka-Volterra" type growth rates (x' — xF(t)) 


Purpose/ 


Invasion: can a mutant invade the 


Evolution of a parametrically 


question 


resident population? 


heterogeneous system over time due to 
natural selection 




Uses both theoretical analysis (bifurcation theory) and numerical solutions 


Visual 


Pairwise invasibility plot (PIP) 


Bifurcation diagram of the corresponding 


representation 




parametrically homogeneous system 



Table 3: A comparison of adaptive dynamics and the Reduction theorem 



Different types of punishment/reward functions 
Case 1. Moderate punishment 

Consider the case, when the punishment function is of the form f(c) = a^pf- This functional form allows 
to incorporate moderate levels of punishment /re ward, the severity of which is determined by the value 
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of the parameter a. 

The pairwise invasibility plot (PIP), which allows to visualize under what relative values of c m and 
c res the dominant eigenvalue A changes its sign, can be seen on Figure [2|3,c, for a = 1 and a = 4. Blue 
regions correspond to the case when A(c res ,c m ) < 0, and consequently the mutant cannot invade; red 
regions correspond to the case when A(c res , c m ) > 0, and the mutant can invade. The point of intersection 
of the two curves corresponds to a convergence stable (CSS) but evolutionarily unstable strategy, which 
can be invaded by "mutants" with large enough values of c m . 

The selection gradient for this punishment function, D(c m ) = 7' — ^ 1+ 2 e a p , so the mutant with c m w c r 
can invade when a<|(l + c m ) 2 (see Figure 0a). 




5 10 5 10 5 10 



a c c 

res res 

Figure 2: Selection gradient and pairwise invasibility plots (PIPs) for function of the type /(c) = Gjt^;- 
Blue regions correspond to parameter values where the mutant cannot invade; red regions correspond to 
parameter values where it can invade, (a) selection gradient, defined in Equation ^ (b) PIP for a = 1; 
the singular strategy c* es is evolutionarily unstable and convergence stable (c) PIP for a = 4; the singular 
strategy c* es is evolutionarily unstable and convergence stable. This punishment/reward function is not 
effective against aggressive over-consumers. 

These results suggest that modest punishment can therefore protect only from modest over-consumers. 
However, more aggressive over-consumers cannot be kept of out the population, as the punishment is not 
severe enough. 

Case 2. Severe punishment /generous reward 

Now consider a case, when the punishment function is of the form /(c) = a(l — c) 3 , which allows to 
impose much more severe punishment on over-consumers and more generous reward on under-consumers 
compared to the previous case. 
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The PIP for this functional form can be seen on Figure [3] for a = 0.1 and a = 0.6. Once again, blue 
regions correspond to the case when A(c res ,c m ) < 0, and consequently the mutant cannot invade; red 
regions correspond to the case when X(c res , c m ) > 0, and the mutant can invade. 

For this type of punishment /reward function, there is a region where invader can invade but unlike 
the previous case, there exists an upper boundary for the possible values of successful c m . Unlike in the 
previous case, singular strategy is stable, which predictably implies that punishment needs to be severe 
enough in order to be able to prevent invasion by over consumers. Moreover, for large enough values of 
a, there can exist two singular strategies, one evolutionarily stable and one evolutionarily unstable. In 
this case, the more "altruistic" ESS, which corresponds to smaller values of c* es , is unstable, probably 
because in this case the reward for under-consumption is insufficiently large, and the punishment is not 
sufficiently severe. The second ESS, which corresponds to larger values of c* es is evolutionarily stable, 
and so the mutant over-consumer cannot invade. 

For this punishment function, the selection gradient is D(c m ) = r — 3a(l — c m ) 2 , so D(c m ) > if 
either c rn > 1 — yg^ or c m < 1 + yg^- One can see that there exists a region, where the mutant cannot 
invade even when the enforcement of punishment, accounted for with parameter a, is quite small. This 
can be interpreted as the rewards of over-consumption being too small below the red invasibility zone, 
and the costs of punishment being too great above the red invasibility zone. 




1 2024 024 



a c c 

res res 

Figure 3: Selection gradient and pairwise invasibility plots for function of the type f(c) — a(l — c) 3 . 
Blue regions correspond to parameter values where the mutant cannot invade; red regions correspond 
to parameter values where it can invade, (a) selection gradient, defined in Equation ^ (b) PIP for 
a = 0.1; the singular strategy c* es is evolutionarily and convergence stable (c) PIP for a = 0.6; there 
seem to appear two singular strategies for large enough values of a; one of the strategies is evolutionarily 
unstable, as neither the reward is large enough for the under-consumers, nor the punishment severe 
enough for over-consumers; the second strategy is c* es is evolutionarily and convergence stable. This 
punishment /reward function is effective against aggressive over-consumers. 
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Case 3. Separating punishment and reward 

Finally, consider the following punishment /reward function: f(c) = p(l — c 71 ). This functional form allows 
to separate the impact of punishment for over-consumption, which is increased or decreased depending 
on the value of parameter 77, and reward for under-consumption, which is influenced primarily by the 
value of parameter p. In this case, the dominant eigenvalue is given by A = r(c m — c res ) — p(c^ — cj? es ). 
As one can clearly see, if for instance r = p, then A > when 77 < 1, and so the mutant should be able to 
invade, since the punishment for over-consumption is always less than its reproductive benefits; if 77 > 1, 
the inequality is reversed, and punishment overwhelms reproductive benefits. If 77 = 1, then everything 
is determined solely by the relative values of growth rate r and the reward parameter p. If c rn w c res , 
then the mutant can (or cannot) invade if D(c m ) = r — prjCm ^ is positive (negative). 

As one can see on Figure [4] and [5] this type of function behaves like case 1 for 77 < 1 and like case 2 
for 77 > 1. These results reiterate the claim that was made in the previous two cases: punishment needs 
to be severe enough in order to successfully prevent invasion by over-consumers. 



punishment/reward function f(c)=p(1-c' 1 ) 




5 10 5 10 5 10 



res res res 

Figure 4: Pairwise invasibility plots for function type /(c) = p(l — c''), r\ = 0.9. The effectiveness of this 
function is the same as was for case 1: punishment is not severe enough to keep over-consumers out of 
the population. 
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Modeling parametrically heterogeneous populations using the Reduction the- 
orem 

Adaptive dynamics allows to obtain analytical results about system susceptibility to invasion in the case 
when the invading "mutant" (in this case an over-consumer with a higher value of c m ) is rare. However, 
what if the initial mutant is not rare (invasion by a group)? What if there is more than one type of 
mutant (the system did not have time to recalibrate before the new mutant came along)? We attempt 
to address these questions through application of the Reduction theorem for replicator equations. 

Without incorporating heterogeneity one cannot study the effects of natural selection on the system, 
and until recently any attempts to write such models resulted in systems of immense dimensionality. 
However, the proposed approach allows us to overcome this problem. 

Assume that each individual in the population studied is characterized by its own value of some intrin- 
sic parameter, such as birth, death, resource consumption rates, etc. (a set of individuals characterized by 
a particular value of the parameter studied is referred to here as a clone) . The distribution of clones will 
change over time due to system dynamics, since different clones impose different selective pressures on 
each other. Consequently, the mean of the parameter studied also changes over time. The mean can be 
computed at each time point using the moment generating function of the initial distribution of clones in 
the population. The changes in the mean of the parameter allows easy tracking of changes in population 
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composition in response to changes in micro-environment (such as varying nutrient availability) or with 
respect to different initial distributions of clones. 

Now let us introduce an auxiliary variable q(t), which satisfies q(t)' = S^y, so that one may rewrite 
the system in the following form: 



g| = r(c -*)<) + /(<=>, 
u s , , s N(t)(l -EHc]) 



Then 



x c {t) = x c (0)e- q ^ +t ^ c+f ^\ (7) 
Total population size is then given by 



N(t) 



f x c (t)dc = 7V(0)e^ (t) f e t(rc+/(c)) P c (0)dc, (8) 
JceA Jc 



where P c (0) = wttti an d the current-time pdf is given by 



c[> N(t) fet(rc+f(c)) Pc ^ dc W 



Now, the expected value of c can be calculated from the definition: 



E l \c\ = / cPJt) = Jc . - —r. ; ' ■ 10 



/ c e*(™+/(c))p c (o)d c ' 

The final system of equations becomes : 



?'(*) = 



p— dz + e 

N(t) 
kz(t) ' 



N(t)(l-Et[c]) 
N(t) + z(t) 



(11) 
(12) 
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where N(t) and -E'[c] are defined above. 

The case, when f(c) = was investigated in (9]. The authors observed that although it takes longer 
for a heterogeneous population to go extinct, tragedy of the commons eventually happens if the maximum 
value of c is large enough. Moreover, one could observe transitional regimes as the population was evolving 
towards being composed of increasingly efficient over-consumers. This situation can be used to forecast 
upcoming crisis and start implementing punishment functions. 

In the proposed form, if the system is parametrically homogeneous system, i.e. when E [c] is constant, 
/(c) can be factored into equation = r (c + /(c) — |f ), and consequently, the spectrum of possible 
dynamical behaviors for this modified system should be qualitatively the same compared to the initial 
model, analyzed in [5j and summarized in Section 2. However, while in [5] nothing prevented increase 
of E [c] up to the maximum possible value, here we want to investigate punishment/reward functions 
that will prevent the increase of E [c] that would otherwise drive the dynamics outside of the regions of 
sustainable coexistence with the common resource. 

Results 

We evaluated the effectiveness of several punishment /reward functions on the dynamics of growth of 
a heterogeneous population of consumers supported by a common dynamical carrying capacity, which 
can be increased or decreased by the individuals themselves. We evaluated the effects of the same 
type of punishment on populations with two different types of initial distributions of clones, namely 
truncated exponential with parameter [i = 10, and Beta distribution with parameters a = 2 , f3 = 2 and 
a = 2, P = 5 (see Figure |6j. Parameter values were chosen in such a way as to give different shapes of the 
initial distribution of clones within the population We hypothesized that the effectiveness of punishment 
will depend not only intrinsic parameter values of the system but also on the initial composition of the 
population. Specifically, we hypothesized that higher levels of punishment /reward will be necessary for 
initial distributions, where population composition is spread out further away from small values of c, such 
as Beta distributions with parameters a = 2, j3 = 5 and even more so Beta distribution with parameters 
a = 2, j3 — 2 (see Figure ^jp) compared to truncated exponential distribution, where fewer over consumers 
are initially present in the population (see Figure pa) . 
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Initial truncated exponential distribution 



Initial Beta distribution 
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Figure 6: Initial distributions (a) truncated exponential with parameter /i = 10 (b) Beta distribution 
with parameters a = 2, f3 = 2 and a = 2, /3 = 5 . 



First, we evaluated the effectiveness of the moderate punishment function of type f(c) = a^pf ; the 
severity of punishment is captured through varying parameter a. The initial distribution was taken to 
be truncated exponential with parameter /i = 10. We took parameter a = 0; 0.5; 1;2 and plotted the 
changes in x c {t) for various c over time (Figure [7]), as well as the changes in the total population size 
and the amount of resource (Figure [8]). We observed that when the punishment imposed is moderate, 
over-consumption could be avoided only when the value of a was very high, i.e. when punishment is 
imposed very severely. 
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200 



Figure 7: (a) Truncated exponential, moderate punishment (case 1). set 1, a = 0. (b) Truncated 
exponential, case 1, set 1, a = 0.5. (c) Truncated exponential, case 1, set 1, a = 1. (d) Truncated 
exponential, case 1, set 1, a = 2. 
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Figure 8: Truncated exponential distribution, set 1, case 1, dynamics of the total population size and total 
resource with respect to different values of a (different levels of severity of imposed punishment). One can 
see that successful management of over consumers was possible only when punishment implementation 
was very high. 



Similar results were observed for the cases of Beta distribution with parameters a — 2, /3 = 2 and 
a = 2, j3 = 5 (additional figures are provided in Appendix). The value of a that was necessary for 
successful management of over-consumers varied depending on different initial distributions, indicating 
that in order to be able to prevent the tragedy of the commons, one needs to evaluate not only the 
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type of punishment and the severity of its enforcement but also match it to the composition of the 
population, since one level of punishment can be effective for one distribution of clones within a population 
of consumers and not another. Noticeably, this kind of insight would be impossible to obtain using the 
analytical methods of adaptive dynamics. 

Next, we conducted the same set of numerical experiments for the severe punishment\generous reward 
function /(c) = a(l — c) 3 . We observed that the value of a that would correspond to successful restraint 
of over-consumers was much lower than in the previous case for all initial distributions considered here 
(see Figures |10[ [9] additional examples are provided in Appendix) . The system was able to support 
individuals with higher values of parameter c present in the initial population than in the previous case. 
In some cases we were also able to observe brief periods of oscillatory transitional dynamical behavior 
before the system collapsed (see Figures [9] and [TO). 
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Figure 9: Severe punishment/generous reward (case 2). Beta distribution with parameters [a, 0] = [2,2]. 
(a) a = 0, (b) a = 0.1, (c) a — 0.17, (d) a — 0.2. One can see the population going through transitional 
regimes before collapse, when the punishment is implemented severely but not quite severely enough 
(when a = 0.7). This most probably corresponds to the expected value of parameter c going through 
region 3 in the phase parameter portrait of the non-distributed system. 
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Figure 10: Case 2, set 1. Beta distribution with parameters [2,2]. Changes in total population size and 
of the common resource over time. 



Finally, we evaluated punishment function of the type /(c) = p(l — c n ), where the intensity of pun- 
ishment and reward are accounted for by parameters r\ and p respectively. We observed that in order to 
evaluate the expected effectiveness of the punishment / reward system one needs to not only adjust param- 



eters p and r\ (see Figure 11 1 to each particular case considered but also be able to evaluate the expected 



range of parameter c (see Figure 12 1 , since one level of punishment may be appropriate for one set of 
initial conditions but not another. For instance, as one can see on Figure [l2j the time to collapse under 
fixed values of parameter p and r/ is different for different initial distributions depending on the maximum 
value of c present in the initial population. Moreover, in accordance with our hypothesis, indeed the 
time to collapse varies depending on the initial distribution of the clones within the population, and the 
higher the frequency of over-consumers is in the initial population, the worse the prognosis. For the 
examples considered, population in which the clones are distributed according to truncated exponential 
distribution is less likely to collapse due to over-consumption than Beta distribution with parameters 
a = 2, f3 = 5 which in turn is slightly less prone to collapse than Beta distribution with parameters 
a = 2, f3 = 2 (Figure |J. 
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Figure 12: The importance of evaluating the range of possible values of c/, illustrated for different initial 
distribution. Punishment function is of the type f(c) — p(l — c v ), where p = 0.6, rj = 1.2. Initial 
distributions are taken to be truncated exponential with parameter p, = 10, and Beta with parameters 
a = 2, j3 = 2 and a = 2, /3 = 5; p — 0.6, r\ = 1.2. The top row corresponds to c e [0, 3]; the bottom row 
corresponds to c € [0,4]. 
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On the existence and properties of stationary distributions 

The adaptive dynamics approach assumes that there exists a stable state of a population composed from 
a single clone x res . In the absence of punishment, any initial distribution evolves in such a way, that 
asymptotically only the clone having maximal possible value of the parameter c will persist. In contrast, 
introducing a punishment function allows for the possibility of a stationary distribution of clones that is 
concentrated in more than one point. Specifically, we show below that a population in a stationary state 
may consist of up to 2 clones with punishment function of the 1st and 3rd type and of up to 3 clones 
with punishment function of the 2nd type. 

Assume that there exists a stationary distribution P c (t) = jf^rj of clones x c that stabilizes over time, 
which would occur when 

= n~^ n ' = Pc{t)[rc + /(c) ~ Et[rc] - Et[m] = °' (13) 

which holds only when 

rc + f(c) = E t [rc]+E t [f(c)]. (14) 



The left hand side of Equation ( 14 1 does not depend on time, while the right hand side is independent 



of c, which implies that Equation (14) holds only when 

rc + f(c)=K, (15) 

where K is a constant. 

Recall that {cj} is the support of a probability distribution P if P(c) > for all Cj and P(c) = 
otherwise. In our case, if there exists a constant K such that 

f(ci) = K-rci, (16) 

then to each K corresponds a stationary distribution P(ci), whose support {cj} coincides with the solution 
to this equation. Therefore, Equation ( |16[ ) can be used to identify the values of c that correspond to 
distributed evolutionary steady states. (Note that a) any distribution concentrated in a single point c 
is clearly stationary, and b) if /(c) = ac, where a is a constant, then every stationary distribution is 
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concentrated in a single point.) 

It follows from Equation (|16[) that for a given punishment function /(c), the set of solutions to Equation 



( 16 1 at fixed K can be represented geometrically as a set S(K) of abscissas of points of intersection of 
the curve y = /(c) and the line y — K — rc. The constant K is free, so by changing K in such a way that 



the set S{K) is not empty, we obtain all possible supports of the stationary distributions (see Figure 13 ) 



V=f(c)=(1-cV(1-tc) 

y=K-rc, K=0.70, r=C4 

y=K-rc, K=0.39, r=0 4 
y=K-rc, K=0.27, r=CA 




— y=f(c)=(1-cy 
- -y=K-rc. K=1.70, r=1 
--y=K-rc, K=0.85. r=1 
y=K-re,K=0.62,M 




Figure 13: Steady states are concentrated at the points of intersection of punishment function y = /(c) 
and a line y = K — rc. Depending on the type of punishment function and on the angle of the line 
y = K — rc, one may have up to three clone types at the stationary state. 



An important conclusion can be made from these figures alone. As one can see in Figure [T3fi , for 
instance, the two points of intersection of /(c) and y = K — rc that correspond to points of concentration 
of clones at a stationary state, can be less than or greater than 1, depending on the parameter r, which 
is the slope of the line y — K — rc. The case when one of the points of intersection is less than one 
("altruistic" under-consuming clones) and one is greater than one (over-consumers) can be interpreted 
as a population of under-consumers "supporting" the population of over-consumers at a stationary state. 
Asymptotic behavior of the population, i.e., population survival or extinction, will be determined by the 



dynamics of System (17 1, which is derived below, and corresponding parameter values. 



Once again, let {c{\ be a support of stationary probability distribution. Then jx Ci = r(K — j^-). 
The right hand side does not depend on i, which implies that 



ln(x Ci (i)) — ln(x Cj (t)) = const, for all i,j 



x (t\ x (0) 

and X H ^ = x Ci (o) for all t. (Note that i, j notation is used in this section to distinguish the fact that 
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we are no longer dealing with rare clones in a resident population, and hence previously used m, res 
notation is no longer appropriate.) 

The dynamics of System Q in this case can be described by 

dN N 

dz eJV . 

where E = E [c] = const and consequently does not depend on t. 

x in) 

Now let Aj = ^"x (o) be the initial frequency of clones x Ci ; then E = J\ A,Cj. Therefore, by changing 
the initial values of x(cj) we can vary the value of E from min(ci) up to max(ci). Stable equilibria (N* , z* ) 



of System (fl7| can now be found from N* — kKz* and 7 — 5z* + » e ^ r * (1 — E) = 7 — Sz* + B = 0, where 



B = ^i+ekf^ ~ cons t- Therefore, z* = given that 7 + B > 0. 



System (17) is very similar to the initial parametrically homogeneous model with no punishment; the 
difference is that now the values of K and E are not identical: the value K is in one-to-one correspondence 
with the support {cj}, and the value of E is defined by the support {c{\ and initial values {x c . (0)}. These 
results suggest that if the system has a stationary distribution, then the dynamics of the resource would 
be determined not by absolute sizes of the populations of 'invaders ' and 'residents ', but by their ratio. 

Remark. In this work we are considering punishment functions, where the support of the stationary 
distribution can consist of up to three points, i.e., the system can be composed of up to three clones (see 
Figure \l3\p ) ■ However, one can conceive of a punishment function that can support an arbitrary number 
of clones in the population at the stationary distribution, which correspond to points of intersection of 
the curve y = /(c) and the line y — K — rc. Such system is of course unlikely to survive for any significant 
period of time unless the number of individuals in the clones with c* 3> 1 is small enough. Such a case can 
be interpreted as a large number of under-consumers supporting a very small number of over-consumers. 



Example 1. 

Consider the punishment function of the form /(c) = "fe^ ; the points of tangency of the curve y — jte^ 
and the line y = K — rc are c\ 2 — ±y // 2a/r — 1. In order for c* = ^Jlajr — 1 > 0, condition 2a > r 
must be satisfied. At c*, f(c*) = yjlar — a, so the line y = K — rc is tangent with y = /(c) = at 
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K = K* = y/2ar — A + r(yj2a/r — 1) (see dotted line in Figure 13 1) 



Hence, the line y = K — rc intersects the curve y = /(c) in two points (see dashed red line in Figure 13 1) 
for all K* <K < a. 



For the purposes of illustration, take a = 1, r = 0.4, K — 0.7 (see Figure 13 1). The stationary 
distribution is concentrated in two points c\ — 0.25 and c 2 = 3, so in its steady state, the population is 
composed both of under-consumers and of over-consumers. The population dynamics is then described 
by system 

dz „ eN , _ x 

at z + N 

where E £ [c\ — 0.25, C2 = 3]. For these particular values of parameters, the system will remain in the 
stable state if Sz* — y + Z *+ N , (1 — E) = 7 + 0.412e(l — E) > 0. Otherwise, it will go to extinction. 

Example 2. 

Now consider the punishment function of the form /(c) = a(l — c) 3 . The line y = K — rc is tangent to 
the curve y = /(c) at c{ 2 = 1 ± y/r/(3a). 

Denote yr/ (3a) = ip. Then c* 2 = 1 ± tp; the condition tp < 1 needs to be satisfied in order for both 
roots to be positive. 

Consequently, 

f(c*, 2) = T?A, 



and therefore 



^1.2 = /(c5 1) + rrf a = r(l ± VO T ^; 



i^min = r(l + V) - ^A 3 < # < r(l - V) + V' 3 = ^ 



There can exist up to 3 points of intersection of the curve /(c) = a(l — c) 3 and the line y = K — rc 
depending on the value of K. 



As one can see in Figure 13 a, where a = r = 1, ip = v/1/3, c^ 2 = 1 ± yl/3, if K > K max = 1.3849, 
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then there are no solutions to Equation ( 16 1; if 1 < K < K max , then there exist 2 solutions to Equation 
(16 1; if 0.6151 = K m i n < K < 1. then there can exist 3 solutions. For instance, taking K — 0.8 and 
solving equation /(c) = K — rc, i.e., (1 — c) 3 = 0.8 — c, we obtain the following solutions: c\ = 0.121, 
C2 = 0.791, C3 = 2.08. So, in this case, at the steady state, two 'altruistic' clones with c\ — 0.121 and 
c-2 = 0.791 are supporting over-consuming clones with C3 = 2.088. 
The overall system dynamics for this example are given by 

f - 

where E e [a = 0.121, c 3 = 2.088]. If k = 1, then at the stationary state N* = 0.8z* and z ,^_" N , = SJ = 
0.444. The system will remain at a steady state if 

eN* 

8z*= 1 + Z . +JV> (1 - £) = 7 + 0.444e(l - > 0; 
Otherwise, it will go to extinction. 



Discussion 

In this paper we studied the dynamics of a consumer-resource type system in order to answer the question 
of whether infliction of punishment for over-consumption and reward for under-consumption can success- 
fully prevent, or at least delay, the onset of the tragedy of the commons. We evaluated the effectiveness 
of three types of punishment functions, as well as the effects of punishment on populations with different 
initial composition of individuals with respect to the levels of resource (over)consumption. 

The proposed model was studied analytically in [51 [TU] without incorporating punishment/reward 
for over-/under- consumption. It describes the interactions of a population of consumers, characterized 
by the value of an intrinsic parameter c, with a common renewable resource in such a way that each 
individual can either contribute to increasing the common dynamical carrying capacity (c < 1) or take 
more than they restore (c > 1). As the value of c increases, the population goes through a series of 
transitional regimes from sustainable coexistence with the resource to oscillatory regime to eventually 
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committing evolutionary suicide through decreasing the common carrying capacity to a level that can no 
longer support the population. 

We began by identifying analytical conditions leading to the possibility of sustainable coexistence 
with the common resource for a subset of cases using adaptive dynamics. This method allows to address 
the question of whether a mutant (in our case, an individual with a higher value of c) can invade a 
homogeneous resident population of resource consumers. We evaluated the effectiveness of three types 
of punishment functions: moderate punishment, /(c) = ayqrf ; severe punishment, /(c) = a(l — c) 3 , 
where the parameter a denotes the severity of implementation of punishment on individuals with the 
corresponding value of parameter c, and function of the type /(c) = p(l — cJ 1 ), which allows to separate 
the influence of reward for under-consumption, primarily accounted for with parameter p. and punishment 
for over-consumption, accounted for with parameter 77. 

We demonstrated that while moderate punishment /reward function can be effective in keeping off 
moderate over-consumers at sufficiently high values of a (see Figures [7] and [8]) , the evolutionarily singular 
strategy in this case is unstable, and thus eventually the population will be invaded by over-consumers 
with large enough value of c. The severe punishment/generous reward approach was uniformly effective, 
almost irrespective of the value of a, and allowed invasion of moderate over consumers only in the small 
region of c w 1, when the punishment is not yet severe enough to outweigh the benefits of moderate over- 
consumption (see Figures [9] and 10 1, and the reward does not yet provide sufficient payoffs in terms of 
higher growth rates. The evolutionarily singular strategy in this case is stable, suggesting that punishment 
is severe enough to prevent the tragedy of the commons. 

Finally, we investigated a punishment/reward function of the type /(c) = p(l — d 1 ) that allowed sep- 
arating the influence of punishment for over-consumption (parameter 77) from that of rewarding under- 
consumption (parameter p). This functional type behaves as a moderate punishment/reward function 
for 77 < 1, rendering it effective only for moderate over-consumers, and it behaves like the severe pun- 
ishment/reward function for 77 > 1; in the critical case 77 = 1 the outcome of the interactions between 
resident and invader populations is determined solely by relative values of growth rate parameter r and 
parameter p. Our results suggest that just rewarding under-consumers is not enough to prevent invasion 
by over consumers and hence one should not expect to be able to prevent the tragedy of the commons 
through reward alone. 

Adaptive dynamics techniques do not yet allow answering the question of system invasibility when 
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the mutant is not rare, such as in cases of migration and consequent invasion by a group, or when 
the resident population is inhomogeneous. We address these questions using the Reduction theorem for 
replicator equations. 

Assume that each individual consumer is characterized by his or her own value of the intrinsic param- 
eter c; a group of consumers with this value of c is referred to here as c-clones. This trait directly affects 
fitness, and thus the distribution of clones will change over time due to system dynamics. The clones 
will experience selective pressures not only from the external environment, competing for the limited 
resource, but also from each other. Consequently, the mean of the parameter will also change over time, 
affecting system dynamics. The mean of the parameter can be computed at each time point from the 
moment generating function of the initial distribution of clones, which allows to evaluate the effective- 
ness of different types of punishment/reward functions on population composition by tracking how the 
distribution of clones changed over time with respect to parameter c. 

We evaluated the effectiveness of the same three types of punishment /reward functions on system 
evolution and calculated predicted the outcomes for different initial distributions of clones within the 
population, which were taken to be truncated exponential with parameter /j = 10, and Beta distribution 
with parameters a — 2, f3 — 2 and a = 2, /3 = 5. The initial distributions were chosen in such a way 
as to give significantly different shapes of the initial probability density function; in applications they 
should be matched to real data, when it is available. We observed that the intensity of implementation 
of punishment\reward has to differ for different initial distributions if one is to successfully stop over- 
consumption, and so in order to be able to make any reasonable predictions one needs to understand 
what the initial composition of the affected population is. We hypothesized that the higher the frequency 
of over consumers in the initial population, the more severe the punishment for over-consumption would 
have to be, and the more generous the reward; specifically, for our examples, we anticipated the prognosis 
to be the most favorable for initial truncated exponential distribution, followed by Beta distribution with 
parameters a = 2, /3 = 5 and then finally a — 2, (3 — 2. This effect additionally implies that the results 
obtained analytically from adaptive dynamics can only be relevant for a subset of cases, i.e. when the 
invader is rare. 

As anticipated, we observed that severe punishment/generous reward approach was much more effec- 
tive in preventing the tragedy of the commons than the moderate punishment/reward function, partic- 
ularly for the cases, when over-consumers were present at higher frequencies (such as both Beta initial 
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distributions). Specifically, we observed that the level of implementation, a, could be nearly ten times 
lower for severe punishment /generous reward system as compared to the moderate punishment function 
(a PS 0.2 vs a ps 2) in order to obtain the same effect of selecting against the over-consumers, which 
can be a very important factor in cases when there are large costs associated with implementation of 
such intervention systems. This comes not only from the severity of punishment but also from the fact 
that moderate punishment allows more time for the over-consumers to replicate, and thus by the time 
the punishment has an appreciable effect, the population composition had changed, and the moderate 
punishment will no longer be effective. So, in punishment implementation one needs to take into account 
not only the severity of punishment but also the time window that moderate punishment may provide, 
allowing over-consumers to proliferate. Within the frameworks of the proposed model, moderate im- 
plementation of more severe punishment\reward system is more effective than severe implementation of 
moderate punishment\reward. 

Proposed here is just one way to try and modify individuals' payoffs in order to prevent resource over- 
consumption - through inflicting punishment and reward that affects the growth rates of clones directly. 
This approach can be modified depending on different situations, inflicting punishment or reward based 
not just on the intrinsic value of c but on total resource currently available. 

Adaptive dynamics and the Reduction theorem. In order to address our questions we have used 
two recently developed methods for modeling parametrically heterogeneous populations: adaptive dynam- 
ics [I] and the Reduction theorem for replicator equations (TJ [5]. Adaptive dynamics allows addressing 
the question of rare "mutant" clone invasion using standard bifurcation theory. The method allows to 
formulate analytical conditions, which can be conveniently visualized using pair wise invasibility plots 
(PIPs). However, the method does not allow addressing questions of system invasibility by clones that 
are not rare, such as in cases of invasion by a group, or when the resident population is heterogeneous. 
These questions can be addressed using the Reduction theorem for replicator equations. 

We consider the moment of "invasion" as the initial time moment; more formally, we assume that all 
"invaders" must be present initially in the population, falling within some known distribution; one can 
then see which clone type(s) will be favored over time due to natural selection and visualize evolutionary 
trajectories, which cannot be achieved using adaptive dynamics. The outcome depends both on the initial 
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distribution of the clones within the population and on the initial state of the system, i.e., other intrinsic 
properties of both individuals within the population and the resources, which in this case were assumed 
to be fixed. However, the system of ODEs that results from the transformation done using the Reduction 
theorem is typically non-autonomous; hence no analytical conditions can typically be obtained using 
standard bifurcation theory. The two methods therefore can complement each other. A more detailed 
comparison of the two methods is in Table [3] 



Tragedy of the commons as prisoner's dilemma 

The conditions that can lead to the tragedy of the commons can be reformulated as a game of prisoner's 
dilemma: in such systems preserving the common resource is in the best interest of the group as a whole 
but over-consumption is in the interest of each particular individual. In the cases where decisions about 
the resource are made by each individual independently, tragedy of the commons seems to be inevitable. 
However, a number of cases have been observed, when the tragedy could be avoided if the individuals 
within the population were able to communicate with each other [14 . The common thread that runs 
through many of the available examples is that 1) in small groups the effects of over-consumption are 
immediately noticeable, since within a small population each individual's actions are more visible than 
in larger populations, and 2) punishment is enforced without delay. 

In this paper we are dealing with a particular case of a game of prisoner's dilemma, where punish- 
ment/reward functions /(c) affect individual payoffs in such a way as to give the possibility to outweigh 



both the benefits of over-consumption and the deleterious effects of under-consumption (see Figure 14). 
We were able to demonstrate that preventing the tragedy of the commons through solely rewarding 
under-consumers is unlikely. However, the effectiveness of "intervention" is increased when both reward 
and punishment systems are enforced, although they do need to be modified depending on population 
composition, since the choices of individuals may be affected not only by their immediate 'payoffs' but 
also by the actions of those surrounding them. Therefore, if one is to try and avert the tragedy of the 
commons through punishing over-consumption and rewarding under-consumption, one needs to not only 
find an appropriate punishment / reward function but also calibrate it to the current composition of each 
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Figure 14: Tragedy of the commons as prisoner's dilemma. The tragedy can be avoided if the immediate 
payoffs of all the players are modified through appropriate punishment/reward functions. 

particular population. 
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